% TABLE9.M

clear variables
close all
clc
rng('default')

oldFolder = pwd;
cd ..
addpath([oldFolder,'/GK18Functions']);
cd(oldFolder);

%% Import data.
% Variables are:
% 1: Real price of gasoline
% 2: CPI Inflation
% 3: 1-yr NSC Inflation expectations
load data;

% Order columns of dataset to be consistent with Definition 3.
% The number of equality restrictions imposed on the ith column of Q (f_i) 
% must satisfy f_1 >= f_2 >= ... >= f_N >= 0. In the case of ties, if the 
% IR of interest is with respect to the jth structural shock, let the jth 
% variable be ordered first.
data = y; clear y; data=data(:,[3 1 2]);

% Set variable names for plotting figures.
%varnames = {'Real Gasoline Price','Inflation','Inflation Expectation'};
varnames = {'Inflation Expectation','Real Gasoline Price','Inflation'};

% Set units for IR. 0 = per cent, 1 = percentage points.
irfUnits = [1 0 0 0 0 0];

%% Options.
opt.p = 12; % No. of lags in VAR
opt.const = 1; % const = 1 if constant in VAR, = 0 otherwise
opt.ivar = [1 2 3];  % Indices of variables whose IR is of interest
opt.cumIR = []; % Indices of variables for which cumulative IR should be used
opt.jshock = 2; % Index of structural shock of interest
opt.H = 16; % Terminal horizon for IR
opt.Kdraws = 1000; % No. of draws from posterior of phi with non-empty IS
opt.aalpha = 0.68; % Credibility level for robust credible and HPD intervals
opt.gridLength = 1000; % Number of points on discrete grid when computing intervals 
opt.signNorm = 1; % = 0 if sign normalisation imposed on A0, = 1 if on A0^(-1)
opt.maxDraw = 3000; % Maximum no. of draws from set of orthonormal matrices
opt.Qdraws = 5; % No. of initial values of Q in optimisation step of Algorithm 1
opt.KTilde = 10000; % No. of draws of Q used in Algorithm 2
opt.dispIter = 100; % Print number of draws remaining every dispIter draws

% Choose algorithm used to compute bounds of IS at each draw from the
% posterior of phi.
% 1 = Algorithm 1 (optimisation) from GK18
% 2 = Algorithm 2 (simulation) from GK18
% 3 = Analytical approach from GMM18
opt.whichAlgo = 2;

% Options for optimisation step in Algorithm 1.
optimOptions = optimoptions('fmincon');
optimOptions.Display = 'off'; % Do not display optimisation output
% Provide analytical gradient of objective function.
optimOptions.GradObj = 'on'; 
% Provide analytical gradient of nonlinear constraint.
optimOptions.SpecifyConstraintGradient = true; 
optimOptions.CheckGradients = false; % Check provided gradients numerically

rng(1); % Set seed for random number generator

%% Input equality and sign restrictions.
% Each row of eqRestr contains a quartet (i,j,l,t) representing a 
% particular equality restriction, where t is the type of restriction:
% t = 1: the (ij)th element of A0 is zero
% t = 2: the (ij)th element of A0^(-1) is zero
% t = 3: the (ij)th element of A_l is zero
% t = 4: the (ij)th element of the long-run cumulative IR (LRCIR) is zero 
% The value of l does not matter except when t = 3.
restr.eqRestr = ...
    [2 1 0 2;  % Real gas price does not respond to idiosyncratic expectations shock on impact
     3 1 0 2]; % Inflation does not respond to idiosyncratic expectations shock on impact

% Each row of signRestr contains a quintet (i,j,h,s,t) representing a
% particular sign restriction, where t is the type of restriction:
% t = 1: the IR of the ith variable with respect to the jth shock at the 
% hth horizon is nonnegative (s = 1) or nonpositive (s = -1).
% t = 2: the (ij)th element of A0 is nonnegative (s = 1) or nonpositive 
% (s = -1). 
% Warning: signRestr should NOT include sign normalisations. These are
% controlled by the 'signNorm' option above. If the sign normalisation is
% placed on A0^(-1), do not place sign restrictions on the response of the
% i-th variable to the i-th shock at h = 0 for i = 1,...,N. Otherwise,
% constraints will be colinear.
%restr.signRestr = ...
%    [2 2 0  1  1;  % Impact response of real gas price to nominal gas price shock is positive
%     3 2 0  1  1;  % Impact response of inflation to nominal gas price shock is positive
%     1 2 0  1  1;  % Impact response of inflation expectation to nominal gas price shock is positive
%     2 3 0 -1  1;  % Impact response of real gas price to core CPI shock is negative
%     3 3 0  1  1;  % Impact response of inflation to core CPI shock is positive
%     1 3 0  1  1;  % Impact response of inflation expectations to core CPI shock is positive
%     1 1 0  1  1];  % Impact response of inflation expectations to idiosyncratic expectations shock is positive
% After dropping diagonal elements:
restr.signRestr = ...
    [3 2 0  1  1;  % Impact response of inflation to nominal gas price shock is positive
     1 2 0  1  1;  % Impact response of inflation expectation to nominal gas price shock is positive
     2 3 0 -1  1;  % Impact response of real gas price to core CPI shock is negative
     1 3 0  1  1];  % Impact response of inflation expectations to core CPI shock is positive


% Conduct posterior inference 
mainfile;
